#include "filters.h"
#include <cstdlib>
#include <ctime>

double greaterModule = 0.0;

double getGreaterModule()
{
    return greaterModule;
}

double getModule(COMPLEX comp)
{
    return sqrt((comp.real*comp.real) + (comp.imag*comp.imag));
}

void periodicNoise( COMPLEX** matrix, int height, int width, int x, int y, float radius, int operation )
{
    float radius_sqr = radius*radius;

    for( int i = 0; i < height; i++ )
    {
        for( int j = 0; j < width; j++ )
        {
            float d1 = (x - j)*(x - j) + (y - i)*(y - i);
            float d2 = (width - x - j)*(width - x - j) + (height - y - i)*(height - y - i);

            if( ( d1 < radius_sqr ) || ( d2 < radius_sqr ) )
            {
                if( operation == PERIODIC_NOISE_OP_REMOVE )
                {
                    matrix[i][j].real = 0;
                    matrix[i][j].imag = 0;
                    matrix[i][j].module = 0;
                }
                else
                {
                    matrix[i][j].real = 180.3125;
                    matrix[i][j].imag = 180.3125;
                    matrix[i][j].module = 255.0;
                }
            }
        }
    }
}

void saltPepperNoise( pixel** matrix, int height, int width, float rate )
{
//    rate /= 100.0;
    srand( time(NULL) );
    for( int i = 0; i < height; i++ )
    {
        for( int j = 0; j < width; j++ )
        {
            double value = (rand()%100);
            //value = value % 100.0;
            if( value < rate )
            {
                double col = rand() % 2;
                matrix[i][j].value = ( col == 0 ) ? 255 : 0;
//                matrix[i][j].value = ( matrix[i][j].value < 128 ) ? 255 : 0;
            }
        }
    }
}

void customNoise( pixel** matrix, int height, int width, float rate )
{
 //   rate /= 100.0;";
    srand( time(NULL) );
    for( int i = 0; i < height; i++ )
    {
        for( int j = 0; j < width; j++ )
        {
            double value = rand() % 100;
//            value /= RAND_MAX;
            if( value < rate )
            {
                double col = rand() % 2;
                int val;
                val = (int)matrix[i][j].value;

                val += ( col == 0 ) ? 40 : -40;
                if(val < 0) val = 0;
                if(val > 255) val = 255;
                matrix[i][j].value = (uByte)val;

            }
        }
    }
}

/** FFT - Fast Fourier Transform */
int Powerof2(int nx, int *m, int *twopm) {
    int pwr;
	*m = 0;
	for (pwr = 1; pwr < nx; pwr = pwr*2) {
		*m = *m + 1;
	}
	*twopm = pwr;
	return true;
}

/*
This function can get a better performance removing the copy operation to and from the real and imag vectors
*/
int FFT2D( COMPLEX** matrixIn, COMPLEX** matrixOut, int nx, int ny, int dir ) {
    int m,twopm;
    double *real,*imag;

    /** Transform the rows */
    real = (double *)malloc(nx * sizeof(double));
    imag = (double *)malloc(nx * sizeof(double));

    if( real == NULL || imag == NULL ) return false;
    if( !Powerof2( nx, &m, &twopm ) || twopm != nx ) return false;

    for( int j = 0; j < ny; j++ )
    {
        for( int i = 0; i < nx; i++ )
        {
            real[i] = matrixIn[i][j].real;
            imag[i] = matrixIn[i][j].imag;
        }
        FFT( dir, m, real, imag );
        for( int i = 0; i < nx; i++ )
        {
            matrixOut[i][j].real = real[i];
            matrixOut[i][j].imag = imag[i];
        }
    }
    free(real);
    free(imag);

    /** Transform the columns */
    real = (double *)malloc(ny * sizeof(double));
    imag = (double *)malloc(ny * sizeof(double));

    if( real == NULL || imag == NULL ) return false ;
    if( !Powerof2( ny, &m, &twopm ) || twopm != ny ) return false;

    for( int i = 0; i < nx; i++ )
    {
        for( int j = 0; j < ny; j++ )
        {
            real[j] = matrixOut[i][j].real;
            imag[j] = matrixOut[i][j].imag;
        }
        FFT( dir, m, real, imag );
        for( int j = 0; j < ny; j++ )
        {
            //double shift = (dir == 1) ? pow(-1.0, i + j) : 1;

            matrixOut[i][j].real = real[j];//*shift;
            matrixOut[i][j].imag = imag[j];//*shift;

            ///Setting the module in each position of the output matrix
            matrixOut[i][j].module = getModule(matrixOut[i][j]);
            if( matrixOut[i][j].module > greaterModule )
            {
                greaterModule = matrixOut[i][j].module;
            }
        }
    }
    free(real);
    free(imag);

    return true;
}

int FFT(int dir,int m,double *x,double *y) {
    long nn,i,i1,j,k,i2,l,l1,l2;
    double c1,c2,tx,ty,t1,t2,u1,u2,z;

    /* Calculate the number of points */
    nn = 1;
    for (i=0;i<m;i++)
        nn *= 2;

    /* Do the bit reversal */
    i2 = nn >> 1;
    j = 0;
    for (i=0;i<nn-1;i++) {
        if (i < j) {
            tx = x[i];
            ty = y[i];
            x[i] = x[j];
            y[i] = y[j];
            x[j] = tx;
            y[j] = ty;
        }
        k = i2;
        while (k <= j) {
            j -= k;
            k >>= 1;
        }
        j += k;
    }

    /* Compute the FFT */
    c1 = -1.0;
    c2 = 0.0;
    l2 = 1;
    for (l=0;l<m;l++) {
        l1 = l2;
        l2 <<= 1;
        u1 = 1.0;
        u2 = 0.0;
        for (j=0;j<l1;j++) {
            for (i=j;i<nn;i+=l2) {
                i1 = i + l1;
                t1 = u1 * x[i1] - u2 * y[i1];
                t2 = u1 * y[i1] + u2 * x[i1];
                x[i1] = x[i] - t1;
                y[i1] = y[i] - t2;
                x[i] += t1;
                y[i] += t2;
            }
            z =  u1 * c1 - u2 * c2;
            u2 = u1 * c2 + u2 * c1;
            u1 = z;
        }
        c2 = sqrt((1.0 - c1) / 2.0);
        if (dir == 1)
            c2 = -c2;
        c1 = sqrt((1.0 + c1) / 2.0);
    }

//    /* Scaling for forward transform */
//    if (dir == 1) {
//        for (i=0;i<nn;i++) {
//            x[i] /= (double)nn;
//            y[i] /= (double)nn;
//        }
//    }
    long n = sqrt(nn);
    for (i=0;i<nn;i++)
    {
        x[i] /= (double)n;
        y[i] /= (double)n;
    }

    return true;
}
